1 Exercice sur la loi Gamma (exercice donné en cours lors de la séance 1)

1.1 Génération d’un échantillon

On génère un échantillon de 50 individus selon une loi Gamma (a,s) où a et s sont tels que : f(x)= 1/( s^a * Gamma(a) ) * x^(a-1) * exp( -(x/s) ) :

On choisit a = 5 et s = 2 :

set.seed(42)
taille <- 50
echantillon <- rgamma(n = taille, shape = 5, scale = 2)

1.2 Maximisation de la fonction maximum de vraissemblance :

Définition de la fonction égale a l’opposée de la vraissemblance :

moins_vraissemblance <- function( p, v ) {
  return(
    - ( 1/( (p[2]^p[1]) * gamma(p[1]) ) ) * ( prod(v) )^( p[1] - 1 ) * exp( -( sum(v) / p[2] ) )
  )
}

Le solveur échoue à trouver les paramètres a et b qui maximisent la vraissemblance. Ici on cherche en fait à minimiser la fonction moins_vraissemblance, ce qui est équivalent. Le solveur renvoie des valeurs différentes selon le point de départ de l’algorithme :

optim(p = c(1,1), fn = moins_vraissemblance, v = echantillon)
## $par
## [1] 1.0 1.1
## 
## $value
## [1] -1.596089e-196
## 
## $counts
## function gradient 
##        3       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

En partant du point (1;1) : a = 1.0 et b = 1.1

optim(p = c(3,1.5), fn = moins_vraissemblance, v = echantillon)
## $par
## [1] 3.0 1.8
## 
## $value
## [1] -2.309985e-25
## 
## $counts
## function gradient 
##        3       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

En partant du point (3;1.5) : a = 3.0 et b = 1.8

optim(p = c(5,4), fn = moins_vraissemblance, v = echantillon)
## $par
## [1]  7.418039 66.840698
## 
## $value
## [1] -1.958337e+288
## 
## $counts
## function gradient 
##      377       NA 
## 
## $convergence
## [1] 10
## 
## $message
## NULL

En partant du point (5;4) : a = 7.42 et b = 66.84

Dans tous les cas, on ne s’approche guère des paramètres réels a = 5 et s = 2. La vraissemblance est difficile à maximiser. Et pourtant, le solveur “optim” fonctionne correctement (sur des cas simples) :

fonction_test <- function( p ) {
  return(
    p[1]^2+p[2]^2
  )
}
optim(p = c(10,10), fn = fonction_test)
## $par
## [1] 0.0003754010 0.0005179101
## 
## $value
## [1] 4.091568e-07
## 
## $counts
## function gradient 
##       63       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Partant du point de coordonnées (10;10), on trouve bien approximativement le minimum de la fonction f(x,y) x²+y², à savoir (0;0).

1.3 Tests asymptotiques

On regarde si l’espérance a x s = 5 x 2 = 10 est dans l’intervalle de confiance asymptotique à 95%, centré en la moyenne empirique, et construit par approximation normale en vertu du théorème central limite. Pour évaluer la pertinence du test, on fait varier n de 2 à 500 et on procède, pour chaque valeur de n, a 100 tests. On observe alors la proportion de test qui accepte l’hypothèse nulle (correcte) : l’espérance vaut 10.

compteur = 0
taux_reussite = 2:500

for( taille in seq(2,500,1) ) {
  
  compteur = 0
  
  for( i in 1:100) {
    
    echantillon <- rgamma(n = taille, shape = 5, scale = 2)
    borne_inf <- mean(echantillon) - ( 1.96 * sd(echantillon) / sqrt(taille) )
    borne_sup <- mean(echantillon) + ( 1.96 * sd(echantillon) / sqrt(taille) )
    
    if( (10 >= borne_inf) & (10 <= borne_sup) ) {
      compteur = compteur + 1
    }
    
    taux_reussite[taille] = compteur
    
  }
}

plot(taux_reussite)

On observe que la fréquence d’acceptation de l’hypothèse nulle (espérance = 10) s’approche très rapidement de la fréquence théorique de 95%, malgré l’approximation normale sur laquelle repose le test.

2 TD0

2.1 Exercice 3 - Loi Weibull

2.1.1 Densité d’une loi de Weibull W()

abscisses <- tibble( abs = seq( 0.01, 10, 0.01 ) )
shape <- tibble( shape = c( 0.5, 1:4 ) )
scale <- tibble( scale = c( 0.5, 1:4 ) )

don <- abscisses %>% crossing(shape) %>% crossing(scale) %>% mutate (ord = dweibull(abs, shape, scale))

ggplot( data = don ) +
  theme_minimal() +
  geom_line( mapping = aes ( x = abs, y = ord) ) +
  facet_grid( shape~scale, scales = "free", labeller = label_both )

2.1.2 Simulation d’une variable aléatoire suivant une loi de Weibull

  • D’abord la méthode simple : on génère un échantillon de 50 valeurs prises par une variable aléatoire suivant la loi de Weibull W(2,1) à l’aide de la fonction déjà implémentée dans R :
echantillon <- rweibull( n= 50, shape = 2, scale = 1)
print(echantillon)
##  [1] 1.17118023 1.18374576 0.59992642 0.88801907 1.09955023 0.49886344
##  [7] 0.90983093 0.98791827 0.16639189 0.33083185 1.18512402 1.53045280
## [13] 1.55964603 0.11662615 0.37898100 0.64306002 1.32625413 0.37208302
## [19] 0.66270373 1.56338206 1.90529765 1.00137051 1.03571506 0.88615928
## [25] 0.84122854 0.87111960 0.70424585 2.22758332 1.16348110 1.21490313
## [31] 1.60532081 1.38245340 0.87019441 1.25212024 2.27621579 1.26418437
## [37] 1.51398186 0.87833235 0.20543892 0.28199209 1.46585914 0.77779746
## [43] 0.09571679 0.90177424 1.22610438 0.88107681 0.92789740 0.82529893
## [49] 1.04263259 0.81960748
  • Ensuite on fait la même chose mais via la méthode du rejet : Ici la densité de notre V.A. suivant la loi W(2,1) est f(x) = 2.x.exp(-x^2). On utilise le fait que f est majorée par c.g(x) où c = 2 et g(x) = exp(-x) (densité d’une loi exponentielle de paramètre 1). Étant donné qu’on sait générer les valeurs de V.A. pour cette loi (et pour la loi uniforme), nous sommes en mesure d’implémenter l’algorithme correspondant à la méthode du rejet :
f <- function(x){
  return( 2*x*exp(-x^2) )
}

g <- function(x){
  return( exp(-x) )
}

#f est majorée par 2g
c <- 2

weibull_sample <- function(n){
  
  sample <- vector( mode = "double", length = n )
  
  for(i in 1:n){
    
    y <- rexp(1)
    u <- runif(1)
    
    while(u > f(y) / ( c * g(y) ) ) {
      y <- rexp(1)
      u <- runif(1)
    }
    
    sample[i] <- y
    
  }
  return (sample)
}

Vérifions que notre algorithme fonctionne en générant 1000 valeurs et en confrontant leur histogramme à la densité théorique :

ggplot( data = tibble( w = weibull_sample(10000) ) ) +
  theme_minimal() +
  geom_histogram( aes( x = w, stat(density)), binwidth = 0.01 ) +
  stat_function( fun = dweibull, args = list( shape = 2,  scale = 1 ), colour = "red", size = 2 )

2.1.3 Convergence de la méthode de simulation

n = 100

valeurs_parametres = c(0.75,1:3)
n_plots = length(valeurs_parametres) * length(valeurs_parametres)
plots = vector("list", n_plots)
i = 1

for(sh in valeurs_parametres){
  for(sc in valeurs_parametres) {
    
    echantillon <- tibble( w = rweibull(100, sh, sc) )
    
    plots[[i]] <- ggplot( data = echantillon ) +
              theme_minimal() +
              ggtitle(str_c("shape = ", sh, " scale = ", sc))+
              geom_histogram( aes( x = w, stat(density)), binwidth = 0.1 ) +
              stat_function( fun = dweibull, args = list( shape = sh,  scale = sc ), colour = "red", size = 1 )
    i = i + 1
  }
}

gridExtra::grid.arrange( grobs = plots, ncol = sqrt(n_plots) )

rm(plots)

n = 1000

plots = vector("list", n_plots)
i = 1

for(sh in valeurs_parametres){
  for(sc in valeurs_parametres) {
    
    echantillon <- tibble( w = rweibull(1000, sh, sc) )
    
    plots[[i]] <- ggplot( data = echantillon ) +
              theme_minimal() +
              ggtitle(str_c("shape = ", sh, " scale = ", sc))+
              geom_histogram( aes( x = w, stat(density)), binwidth = 0.1 ) +
              stat_function( fun = dweibull, args = list( shape = sh,  scale = sc ), colour = "red", size = 1 )
    i = i + 1
  }
}

gridExtra::grid.arrange( grobs = plots, ncol = sqrt(n_plots) )

rm(plots)

n = 10 000

plots = vector("list", n_plots)
i = 1

for(sh in valeurs_parametres){
  for(sc in valeurs_parametres) {
    
    echantillon <- tibble( w = rweibull(10000, sh, sc) )
    
    plots[[i]] <- ggplot( data = echantillon ) +
              theme_minimal() +
              ggtitle(str_c("shape = ", sh, " scale = ", sc))+
              geom_histogram( aes( x = w, stat(density)), binwidth = 0.1 ) +
              stat_function( fun = dweibull, args = list( shape = sh,  scale = sc ), colour = "red", size = 1 )
    i = i + 1
  }
}

gridExtra::grid.arrange( grobs = plots, ncol = sqrt(n_plots) )

rm(plots)

La convergence semble se faire d’autant plus lentement que le premier paramètre (shape) s’approche de 0.

2.1.4 Anecdote sur la loi de Weibull

La loi de Weibull est utilisée pour évaluer la durée de vie de systèmes dont la probabilité de défaillir est une fonction puissance du temps.

3 TD0

3.1 Exercice 4 - Loi des Grands Nombres

3.1.1 Étude de la loi uniforme

On génère 1000 échantillons de n variables aléatoires X suivant une loi uniforme. Puis on représente l’histogramme de la V.A. sqrt(n) * ( Sn/n - E(X) ) / sigma(X), ainsi que la densité de la loi normale

v <- 1:1000
taille <- c(10, 100, 1000, 10000)

for (n in taille) {
  for (i in 1:1000) {
    v[i]  <- sqrt(12*n) * ((sum(runif(n)) / n) - 0.5)
  }
  
  echantillon <- tibble(unif = v)
  
  plot <- ggplot(data = echantillon) +
    theme_minimal() +
    ggtitle(str_c("n = ", n)) +
    geom_histogram(aes(x = unif, stat(density)), binwidth = 0.05) +
    stat_function(fun = dnorm,
                  colour = "red",
                  size = 1)
  
  print(plot)
}

On observe la convergence vers l’espérance (Loi des Grands Nombres) de la moyenne empirique ainsi que son comportement normal autour de l’espérance.

3.1.2 Même chose pour la loi de Poisson de paramètre 1 :

for (n in taille) {
  for (i in 1:1000) {
    v[i]  <- sqrt(n) * ((sum(rpois(n, lambda = 1)) / n) - 1)
  }
  
  echantillon <- tibble(pois = v)
  
  plot <- ggplot(data = echantillon) +
    theme_minimal() +
    ggtitle(str_c("n = ", n)) +
    geom_histogram(aes(x = pois, stat(density)), binwidth = 0.05) +
    stat_function(fun = dnorm,
                  colour = "red",
                  size = 1)
  
  print(plot)
}

4 Cartes de contrôle

4.1 Carte X, R

On génère deux échantillons E1 et E2 de grande taille (100 000) :
- le premier échantillon suit une loi normale d’espérance 1 / écart-type 0.5
- le deuxième échantillon suit une loi normale d’espérance 1.1 / écart-type 1

On peut imaginer que cette situation correspond à une machine qui se dérègle de la façon suivante : la moyenne de la variable change un peu et la stabilité de la production diminue beaucoup au sens où l’écart-type augmente sensiblement.

set.seed(1)
sample1 <- rnorm(100000, 1, 0.5)
sample2 <- rnorm(100000, 1.1, 1)

Supposons que l’on ne connaisse aucun paramètre. Tout d’abord testons la normalité de l’échantillon, via un sous-échantillon de taille 30 :

shapiro.test(sample(sample1, 30))
## 
##  Shapiro-Wilk normality test
## 
## data:  sample(sample1, 30)
## W = 0.9911, p-value = 0.9955

On accepte la normalité de l’échantillon. Nous sommes bien dans le cadre théorique des cartes de contrôle. On estime la moyenne X et l’étendue R à partir de 25 échantillons de taille 10. Ces données joueront le rôle de consigne dans les cartes de contrôle. On calcule la moyenne et l’étendue de chacun des 25 échantillons, et on prend la moyenne de ces 25 valeurs pour chaque paramètre.

x_estimations <- vector("double", 25)
R_estimations <- vector("double", 25)
for(i in 1:25)
{
  little_sample <- sample(sample1, 10)
  x_estimations[i] <- mean(little_sample)
  R_estimations[i] <- max(little_sample) - min(little_sample)
}

x <- mean(x_estimations)
R <- mean(R_estimations)

On trouve les coefficient A2, D3 et D4 et dans les tables pour des échantillons de taille 10. On en déduit les LCL/UCL. On prend les estimations X et R comme consigne.

a2 <- 0.308
d3 <- 0.223
d4 <- 1.777
uclx <- x + a2 * R
lclx <- x - a2 * R
uclr <- d4 * R
lclr <- d3 * R

À présent on simule le suivi de la moyenne et de l’étendue sur 30 échantillons, tous de taille 10.
- Les 15 premiers échantillons proviennent du super-échantillon E1 (loi normale d’espérance 1 / écart-type 0.5)
- Les 15 échantillons suivants proviennent du super-échantillon E2 (loi normale d’espérance 1.1 / écart-type 0.5)

Cette situation peut correspondre à un problème sur une ligne de production, un dérèglement de machine, une erreur humaine, etc…

nombre_observation <- 30
x_samples <- vector("double", nombre_observation)
r_samples <- vector("double", nombre_observation)
for(i in 1:30)
{
  if(i <= 15)
  {
    sample <- sample(sample1, 10)
    x_samples[i] <- mean(sample)
    r_samples[i] <- max(sample) - min(sample)
  }else
  {
    sample <- sample(sample2, 10)
    x_samples[i] <- mean(sample)
    r_samples[i] <- max(sample) - min(sample)
  }
}

plotx <- ggplot(data = tibble(moyenne = x_samples)) +
  aes( x = 1:nombre_observation, y = moyenne) +
  geom_line() +
  geom_hline(yintercept = x, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclx, color = "red", size = 1) +
  geom_hline(yintercept = uclx, color = "red", size = 1)+
  geom_text(x = 28, y = x + 0.05, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 28, y = lclx + 0.05, label = "LCL", color = "red") +
  geom_text(x = 28, y = uclx - 0.05, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("Moyenne") +
  theme_minimal() +
  ggtitle("Contrôle de la moyenne X\u0305") +
  theme( plot.title = element_text(hjust = 0.5) )

plotr <- ggplot(data = tibble(etendue = r_samples)) +
  aes( x = 1:nombre_observation, y = etendue) +
  geom_line() +
  geom_hline(yintercept = R, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclr, color = "red", size = 1) +
  geom_hline(yintercept = uclr, color = "red", size = 1)+
  geom_text(x = 6, y = R + 0.25, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 6, y = lclr + 0.25, label = "LCL", color = "red") +
  geom_text(x = 6, y = uclr + 0.25, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("Étendue") +
  theme_minimal() +
  ggtitle("Contrôle de l'étendue R") +
  theme( plot.title = element_text(hjust = 0.5) )
  
print(plotx)

print(plotr)

On constate que la moyenne a probablement augmenté même si la carte de contrôle correspondante ne produit pas d’alerte. En revanche, l’étendue n’est plus du tout sous contrôle.

4.2 Carte X, s

On génère deux échantillons E1 et E2 de grande taille (100 000) :
- le premier échantillon suit une loi normale d’espérance 1 / variance 1
- le deuxième échantillon suit une loi normale d’espérance 1.5 / variance 0.5

On peut imaginer que cette situation correspond à une machine qui se dérègle de la façon suivante : la moyenne de la variable est modifiée et prend une valeur non-conforme et la stabilité de la production augmente au sens où la variance diminue sensiblement. La machine se met à produire beaucoup de pièces non-conformes.

set.seed(42)
sample1 <- rnorm(100000, 1, 1)
sample2 <- rnorm(100000, 1.8, 0.8)

Supposons que l’on ne connaisse aucun paramètre. Tout d’abord testons la normalité de l’échantillon, via un sous-échantillon de taille 25 :

shapiro.test(sample(sample1, 25))
## 
##  Shapiro-Wilk normality test
## 
## data:  sample(sample1, 25)
## W = 0.95142, p-value = 0.2699

On accepte la normalité de l’échantillon. Nous sommes bien dans le cadre théorique des cartes de contrôle. On estime la moyenne X et l’écart-type s à partir de 10 échantillons de taille 25. En effet la carte de contrôle X, s est indiquée lorsque la taille des sous-échantillons est supérieure à 10 ou 12, car lorsque la taille augmente l’étendue perd en efficacité. Ces données joueront le rôle de consigne dans les cartes de contrôle. On calcule la moyenne et l’écart-type de chacun des 10 échantillons, et on prend la moyenne de ces 10 valeurs pour chaque paramètre.

Remarque : la fonction sd() de R utilise bien le dénominateur n-1.

x_estimations <- vector("double", 10)
s_estimations <- vector("double", 10)
for(i in 1:10)
{
  little_sample <- sample(sample1, 25)
  x_estimations[i] <- mean(little_sample)
  s_estimations[i] <- sd(little_sample)
}

x <- mean(x_estimations)
s <- mean(s_estimations)

On trouve les coefficient A3, B3 et B4 et dans les tables pour des échantillons de taille 25. On en déduit les LCL/UCL. On prend les estimations X et s comme consigne.

a3 <- 0.606
b3 <- 0.565
b4 <- 1.435
uclx <- x + a3 * s
lclx <- x - a3 * s
ucls <- b4 * s
lcls <- b3 * s

À présent on simule le suivi de la moyenne et de l’étendue sur 30 échantillons, tous de taille 25
- Les 15 premiers échantillons proviennent du super-échantillon E1 (loi normale d’espérance 1 / variance 1)
- Les 15 échantillons suivants proviennent du super-échantillon E2 (loi normale d’espérance 1.8 / variance 0.8)

Cette situation peut correspondre à un problème sur une ligne de production, un dérèglement de machine, une erreur humaine, etc…

nombre_observation <- 30
x_samples <- vector("double", nombre_observation)
s_samples <- vector("double", nombre_observation)
for(i in 1:30)
{
  if(i <= 15)
  {
    sample <- sample(sample1, 25)
    x_samples[i] <- mean(sample)
    s_samples[i] <- sd(sample)
  }else
  {
    sample <- sample(sample2, 25)
    x_samples[i] <- mean(sample)
    s_samples[i] <- sd(sample)
  }
}

plotx <- ggplot(data = tibble(moyenne = x_samples)) +
  aes( x = 1:nombre_observation, y = moyenne) +
  geom_line() +
  geom_hline(yintercept = x, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclx, color = "red", size = 1) +
  geom_hline(yintercept = uclx, color = "red", size = 1)+
  geom_text(x = 25, y = x + 0.1, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 25, y = lclx + 0.1, label = "LCL", color = "red") +
  geom_text(x = 25, y = uclx + 0.1, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("Moyenne") +
  theme_minimal() +
  ggtitle("Contrôle de la moyenne X\u0305") +
  theme( plot.title = element_text(hjust = 0.5) )

plots <- ggplot(data = tibble(sd = s_samples)) +
  aes( x = 1:nombre_observation, y = sd) +
  geom_line() +
  geom_hline(yintercept = s, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lcls, color = "red", size = 1) +
  geom_hline(yintercept = ucls, color = "red", size = 1)+
  geom_text(x = 25, y = s + 0.05, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 25, y = lcls + 0.05, label = "LCL", color = "red") +
  geom_text(x = 25, y = ucls - 0.05, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("s") +
  theme_minimal() +
  ggtitle("Contrôle de l'écart type s") +
  theme( plot.title = element_text(hjust = 0.5) )
  
print(plotx)

print(plots)

On constate que l’écart-type est encore sous contrôle malgré une tendance à la baisse visible. En revanche, la moyenne n’est plus du tout sous contrôle.

4.3 Carte p

On génère un échantillon de grande taille (1 000 000) suivant une loi de Bernoulli de paramètre 0.03. Cette situation pourrait être celle d’une machine qui fabrique des pièces en grande quantité, telle que 3 pièces sur 100 sont défectueuses. On génère un autre échantillon de même taille, mais où suite à un problème donné, le taux de non-conformité de la machine est passé à 10 pièces sur 100 :

set.seed(42)
big_sample1 <- rbinom(n = 1000000, size = 1, prob = 0.03)
big_sample2 <- rbinom(n = 1000000, size = 1, prob = 0.10)

Supposons que l’on ne connaisse pas la probabilité de défaut. On estime la moyenne p à partir de m = 30 sous-échantillons de taille n = 100. Les sous-échantillons proviennent de l’échantillon big_sample1, dans lequel on estime que le taux de non-conformité est sous contrôle. Cette donnée joue le rôle de consigne dans la carte de contrôle. On calcule la moyenne des taux de non-conformité associés à chacun des sous-échantillons, ce qui revient dans ce cas précis (p chart) à calculer directement le taux de non-conformité sur les n x m pièces :

m <- 30
n <- 100
p <- mean( sample( big_sample1, m*n ) )

On en déduit les LCL/UCL. On prend l’estimation de p comme consigne.

if(p - 3 * sqrt( p * (1-p) / n ) < 0){
  
}
## NULL
lclp <- p - 3 * sqrt( p * (1-p) / n )
uclp <- p + 3 * sqrt( p * (1-p) / n )

À présent on simule le suivi du taux de non-conformité sur 30 sous-échantillons, tous de taille n = 100.
- Les 15 premiers échantillons proviennent de l’échantillon big_sample1 (taux de non-conformité réel : 0.03)
- Les 15 échantillons suivants proviennent de l’échantillon big_sample2 (taux de non-conformité réel : 0.10)

Cette situation peut correspondre à un problème sur une ligne de production, un dérèglement de machine, une erreur humaine, etc…

nombre_observation <- 30
p_samples <- vector("double", nombre_observation)
for(i in 1:30)
{
  if(i <= 15)
  {
    p_samples[i] <- mean(sample(big_sample1, n))
  }else
  {
    p_samples[i] <- mean(sample(big_sample2, n))
  }
}

plotp <- ggplot(data = tibble(moyenne = p_samples)) +
  aes( x = 1:nombre_observation, y = moyenne) +
  geom_line() +
  geom_hline(yintercept = p, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclp, color = "red", size = 1) +
  geom_hline(yintercept = uclp, color = "red", size = 1)+
  geom_text(x = 25, y = p + 0.01, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 25, y = lclp + 0.01, label = "LCL", color = "red") +
  geom_text(x = 25, y = uclp + 0.01, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("p") +
  theme_minimal() +
  ggtitle("Contrôle du taux de non-conformité p") +
  theme( plot.title = element_text(hjust = 0.5) )
  
print(plotp)

On constate que le taux de non-conformité n’est plus sous contrôle à partir du sous-échantillon n°16.

4.4 Carte c

Supposons que le nombre moyen de pièces défectueuses par sous-échantillon suive une loi de Poisson de paramètre connu c = 15. Cette donnée joue le rôle de consigne dans la carte de contrôle. On calcule donc directement les LCL/UCL :

c <- 15
lclc <- c - 3 * sqrt( c )
uclc <- c + 3 * sqrt( c )

À présent on simule le nombre de pièces défectueuses sur 30 sous-échantillons.
- Les 15 premiers sous-échantillons présentent 15 pièces défectueuses en moyenne
- Les 15 échantillons suivants présentent 32 pièces défectueuses en moyenne

Cette situation peut correspondre à un problème sur une ligne de production, un dérèglement de machine, une erreur humaine, etc…

nombre_observation <- 30
c_data <- tibble( c_data = c( rpois(15,15), rpois(15,32) ) )

plotc<- ggplot(data = c_data) +
  aes( x = 1:nombre_observation, y = c_data) +
  geom_line() +
  geom_hline(yintercept = c, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclc, color = "red", size = 1) +
  geom_hline(yintercept = uclc, color = "red", size = 1)+
  geom_text(x = 27, y = c + 2, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 27, y = lclc + 2, label = "LCL", color = "red") +
  geom_text(x = 27, y = uclc + 2, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("c") +
  theme_minimal() +
  ggtitle("Contrôle du nombre de pièces défectueuses") +
  theme( plot.title = element_text(hjust = 0.5) )
  
print(plotc)

Le nombre de pièces défectueuses n’est plus sous contrôle à partir du sous-échantillon n°16 environ.

4.5 Carte g,h

Supposons que la proportion moyenne pièces défectueuses soit p = 0.03. On suppose ce paramètre inconnu lors de la réalisation de la carte de contrôle. Une opération de contrôle consiste ici à inspecter des pièces jusqu’à trouver n = 50 pièces défectueuses. On suppose que le nombre de pièces à inspecter avant de tomber sur une non-conformité suit une loi géométrique de paramètre p = 0.03.

Tout d’abord estimons non pas p mais tn, le nombre moyen de pièces à inspecter avant de trouver n pièces défectueuses. Pour cela on utilise l’estimateur tn, la moyenne des tni (nombre de pièces inspectées jusqu’à trouver n pièces défectueuses) pour 0 ≤ i ≤ m, pù m = 50.

p <- 0.03
n <- 50
m <- 50
t_estimations <- vector("double", m)
for(i in 1:m)
{
  number_events <- sum( rgeom( n, p ) )
  t_estimations[i] <- number_events
}
t <- mean(t_estimations)

On en déduit les valeurs des cartes de contrôle g et h :

clg <- t
uclg <- t + 3 * sqrt( n * ( (t/n) - 1 ) * (t/n) )
lclg <- t - 3 * sqrt( n * ( (t/n) - 1 ) * (t/n) )
  
clh <- t/n
uclh <- (t/n) + (3/sqrt(n)) * sqrt( ( (t/n) - 1 ) * (t/n) )
lclh <- (t/n) - (3/sqrt(n)) * sqrt( ( (t/n) - 1 ) * (t/n) )

À présent on simule 30 valeurs de tn (nombre de pièces inspectées avant de trouver n non-conformités) :
- Pour les 15 premières valeurs : p = 0.03 comme ci-dessus
- Pour les 15 valeurs suivantes : p’ = 0.08 (augmentation soudaine du taux de non-conformité suite à un problème donné sur la ligne de production)

nombre_observation <- 30
t_values <- vector("double", nombre_observation)

for(i in 1:30)
{
  if(i <= 15)
  {
    t_values[i] <- sum( rgeom( n, p ) )
  }else
  {
    t_values[i] <- sum( rgeom( n, 0.08 ) )
  }
}

plotg <- ggplot(data = tibble(t_data = t_values) ) +
  aes( x = 1:nombre_observation, y = t_data) +
  geom_line() +
  geom_hline(yintercept = clg, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclg, color = "red", size = 1) +
  geom_hline(yintercept = uclg, color = "red", size = 1)+
  geom_text(x = 25, y = clg + 50, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 25, y = lclg + 50, label = "LCL", color = "red") +
  geom_text(x = 25, y = uclg + 50, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("t") +
  theme_minimal() +
  ggtitle("g chart") +
  theme( plot.title = element_text(hjust = 0.5) )

ploth <- ggplot(data = tibble(t_data = t_values/n) ) +
  aes( x = 1:nombre_observation, y = t_data) +
  geom_line() +
  geom_hline(yintercept = clh, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclh, color = "red", size = 1) +
  geom_hline(yintercept = uclh, color = "red", size = 1)+
  geom_text(x = 25, y = clh + 50/n, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 25, y = lclh + 50/n, label = "LCL", color = "red") +
  geom_text(x = 25, y = uclh + 50/n, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("t/n") +
  theme_minimal() +
  ggtitle("h chart") +
  theme( plot.title = element_text(hjust = 0.5) )


  
print(plotg)

print(ploth)

Le nombre de pièces défectueuses n’est effectivement plus sous contrôle à partir du sous-échantillon n°16.

4.6 CUSUM

4.6.1 CUSUM avec phase préliminaire d’estimation des paramètres / Suivi de la production via une carte de contrôle de Shewhart

On génère deux échantillons E1 et E2 de grande taille (100 000) :
- le premier échantillon suit une loi normale d’espérance 5 / variance 1
- le deuxième échantillon suit une loi normale d’espérance 5.3 / variance 1

On peut imaginer que cette situation correspond à une machine qui produit des pièces de taille 5 en moyenne, avec un écart-type 1. Suite à un problème donné, la machine commence à produire des pièces de taille 5.3 en moyenne, avec le même écart-type 1. Dans un premier temps, on se propose de suivre l’évolution de la production avec une carte de contrôle de Shewhart.

set.seed(4)
big_sample1 <- rnorm(100000, 5, 1)
big_sample2 <- rnorm(100000, 5.3, 1)

Supposons que l’on ne connaisse aucun paramètre. Tout d’abord testons la normalité de l’échantillon, via un sous-échantillon de taille 25 :

shapiro.test(sample(big_sample1, 30))
## 
##  Shapiro-Wilk normality test
## 
## data:  sample(big_sample1, 30)
## W = 0.96583, p-value = 0.4321

On accepte la normalité de l’échantillon. Nous sommes bien dans le cadre théorique de la carte de contrôle de Shewhart. On estime la moyenne X et l’écart-type s à partir de 50 échantillons de taille 10. Ces données joueront le rôle de consigne dans la carte de contrôle. On calcule la moyenne et l’écart-type de chacun des 50 échantillons, et on prend la moyenne de ces 50 valeurs pour chaque paramètre.

x_estimations <- vector("double", 50)
s_estimations <- vector("double", 50)
for(i in 1:50)
{
  little_sample <- sample(big_sample1, 10)
  x_estimations[i] <- mean(little_sample)
  s_estimations[i] <- sd(little_sample)
}

x <- mean(x_estimations)
s <- mean(s_estimations)

On trouve les coefficient A3, B3 et B4 et dans les tables pour des échantillons de taille 10. On en déduit les LCL/UCL. On prend les estimations X et s comme consigne.

a3 <- 0.975
b3 <- 0.284
b4 <- 1.716
uclx <- x + a3 * s
lclx <- x - a3 * s
ucls <- b4 * s
lcls <- b3 * s

À présent on simule le suivi de la moyenne et de l’étendue sur 30 échantillons, tous de taille 10.
- Les 15 premiers échantillons proviennent du super-échantillon E1 (loi normale d’espérance 5 / variance 1)
- Les 15 échantillons suivants proviennent du super-échantillon E2 (loi normale d’espérance 5.3 / variance 1)

nombre_observation <- 30
x_samples <- vector("double", nombre_observation)
s_samples <- vector("double", nombre_observation)
for(i in 1:30)
{
  if(i <= 15)
  {
    sample <- sample(big_sample1, 10)
    x_samples[i] <- mean(sample)
    s_samples[i] <- sd(sample)
  }else
  {
    sample <- sample(big_sample2, 10)
    x_samples[i] <- mean(sample)
    s_samples[i] <- sd(sample)
  }
}

plotx <- ggplot(data = tibble(moyenne = x_samples)) +
  aes( x = 1:nombre_observation, y = moyenne) +
  geom_line() +
  geom_hline(yintercept = x, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclx, color = "red", size = 1) +
  geom_hline(yintercept = uclx, color = "red", size = 1)+
  geom_text(x = 1, y = x + 0.1, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 1, y = lclx + 0.1, label = "LCL", color = "red") +
  geom_text(x = 1, y = uclx - 0.1, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("Moyenne") +
  theme_minimal() +
  ggtitle("Contrôle de la moyenne X\u0305") +
  theme( plot.title = element_text(hjust = 0.5) )

plots <- ggplot(data = tibble(sd = s_samples)) +
  aes( x = 1:nombre_observation, y = sd) +
  geom_line() +
  geom_hline(yintercept = s, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lcls, color = "red", size = 1) +
  geom_hline(yintercept = ucls, color = "red", size = 1)+
  geom_text(x = 1, y = s + 0.05, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 1, y = lcls + 0.05, label = "LCL", color = "red") +
  geom_text(x = 1, y = ucls - 0.05, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("s") +
  theme_minimal() +
  ggtitle("Contrôle de l'écart type s") +
  theme( plot.title = element_text(hjust = 0.5) )
  
print(plotx)

print(plots)

La carte de contrôle de Shewhart ne détecte aucune anomalie. Au vu des graphiques, on suspecte néanmoins une hausse de la moyenne X. On se propose donc de réaliser une carte de contrôle CUSUM pour confirmer ou infirmer cette hypothèse. Les paramètres réels étant inconnus, on affine nos estimations avec les échantillons de la carte de contrôle (du moins ceux où la production semble être sous contrôle) :

x_new <- mean(c(x_estimations, x_samples[1:15]))
s_new <- mean(c(s_estimations, s_samples[1:15]))

Notre estimation de X passe de 5.02 à x = 5 et celle de s de 0.94 à s = 0.95. Supposons que le décalage de X que nous souhaitons détecter soit delta = 0.2. On pose K = delta . s / 2 et H = 5.s.

x <- x_new
s <- s_new
d <- 0.25
k <- d*s/2
h = 5*s

On génère 30 pièces : 12 pièces de l’échantillon 1 (loi N(5,1)), 18 pièces de l’échantillon 2 (loi N(5.3,1)). On construit la carte CUSUM associée.

nbre <- 30
x_values <- c(sample(big_sample1,12), sample(big_sample2,18))
cp <- vector("double", nbre)
cm <- vector("double", nbre)
np <- vector("double", nbre)
nm <- vector("double", nbre)

cp[1] <- max(0, x_values[1] - (x+k) )
cm[1] <- max(0, (x+k) - x_values[1] )
np[1] <- if_else(cp[1]==0, 0, 1)
nm[1] <- if_else(cm[1]==0, 0, 1)
  
for(i in 2:30){
  cp[i] <- max(0, x_values[i] - (x+k) + cp[i-1] )
  cm[i] <- max(0, (x+k) - x_values[i] + cm[i-1] )
  np[i] <- if_else(cp[i]==0, 0, np[i-1] + 1 )
  nm[i] <- if_else(cm[i]==0, 0, nm[i-1] + 1 )
}

cusum <- tibble(i = 1:nbre, x_i = round(x_values,2), C_plus = round(cp,2), N_plus = np, C_moins = round(cm,2), N_moins = nm)

kable(cusum)
i x_i C_plus N_plus C_moins N_moins
1 5.15 0.03 1 0.00 0
2 5.24 0.15 2 0.00 0
3 3.65 0.00 0 1.47 1
4 4.48 0.00 0 2.11 2
5 6.00 0.88 1 1.22 3
6 4.82 0.59 2 1.52 4
7 3.67 0.00 0 2.97 5
8 6.04 0.92 1 2.05 6
9 4.13 0.00 0 3.04 7
10 5.64 0.53 1 2.51 8
11 6.53 1.94 2 1.10 9
12 4.38 1.20 3 1.84 10
13 4.97 1.06 4 1.98 11
14 6.73 2.67 5 0.37 12
15 5.35 2.90 6 0.14 13
16 5.60 3.38 7 0.00 0
17 5.90 4.17 8 0.00 0
18 4.55 3.60 9 0.57 1
19 6.67 5.15 10 0.00 0
20 4.95 4.98 11 0.17 1
21 4.40 4.26 12 0.89 2
22 5.16 4.30 13 0.85 3
23 5.63 4.81 14 0.34 4
24 6.39 6.08 15 0.00 0
25 6.55 7.51 16 0.00 0
26 6.97 9.37 17 0.00 0
27 5.50 9.74 18 0.00 0
28 5.15 9.77 19 0.00 0
29 5.99 10.65 20 0.00 0
30 5.97 11.50 21 0.00 0 
plot_cusum <- ggplot( data = cusum %>% select(i, C_plus, C_moins) %>% mutate(C_moins = - C_moins) %>% gather(type, value, C_plus:C_moins) ) +
  aes( x = i, y = value, color = type) +
  geom_line(size = 1) +
  geom_hline(yintercept = 0, color = "black", size = 1) +
  geom_hline(yintercept = h, color = "red", size = 1.2) +
  geom_hline(yintercept = -h, color = "red", size = 1.2)+
  geom_text(x = 1, y = h + 1, label = "H", color = "red") +
  geom_text(x = 1, y = -h + 1, label = "-H", color = "red") +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:30 ) +
  ylab("C+ / C-") +
  theme_minimal() +
  ggtitle("CUSUM chart") +
  theme( plot.title = element_text(hjust = 0.5) )

print(plot_cusum)

On constate sur le graphique qu’on dépasse la valeur de H pour la pièce i = 19. En se référant au tableau CUSUM, on lit N+ = 10 lorsque i = 19. On peut donc situer l’évènement (erreur humaine, dérèglement ou autre) qui a provoqué le décalage de la moyenne aux alentours de la pièce i = 9. Cette inférence est satisfaisante, quand on compare à l’origine réelle du problème en i = 12.

4.6.2 CUSUM autonome

On travaille sur deux échantillons :
- le premier échantillon suit une loi normale d’espérance 5 / variance 1
- le deuxième échantillon suit une loi normale d’espérance 5.5 / variance 1

On génère 12 observations issues du premier échantillon et 18 observations issues du second. On construit directement la carte CUSUM autonome proposée par Hawkins.

nbre <- 30
set.seed(42)
big_sample1 <- rnorm(100000, 5, 1)
big_sample2 <- rnorm(100000, 5.5, 1)
x_values <- c(sample(big_sample1,12), sample(big_sample2,18))

x <- vector("double", nbre)
w <- vector("double", nbre)
s <- vector("double", nbre)
t <- vector("double", nbre)
u <- vector("double", nbre)

x[1] <- x_values[1]
w[1] <- 0
s[1] <- 0
t[1:2] <- 0
u[1:2] <- 0
  
for(i in 2:30){
  x[i] <- x[i-1] + ( ( x_values[i] - x[i-1] ) / i )
  w[i] <- w[i-1] + (  (i-1) * (x_values[i] - x[i-1])*(x_values[i] - x[i-1]) / i  )
  s[i] <- sqrt( w[i] / (i-1) )
}

for(i in 3:30){
  t[i] <- ( x_values[i] - x[i-1] ) / s[i-1]
  a <- sqrt( (i-1) / i )
  u[i] <- qnorm( pt( a*t[i], i-2 ) )
}

# CUSUM sur les valeurs de U ~ N(0,1)

u <- u[3:nbre]
k <- 0.1
h <- 5
mu <- 0
sigma <- 1
cp <- vector("double", nbre-2)
cm <- vector("double", nbre-2)
np <- vector("double", nbre-2)
nm <- vector("double", nbre-2)

cp[1] <- max(0, u[1] - (mu+k) )
cm[1] <- max(0, (mu+k) - u[1] )
np[1] <- if_else(cp[1]==0, 0, 1)
nm[1] <- if_else(cm[1]==0, 0, 1)
  
for( i in 2:(nbre-2) ){
  cp[i] <- max(0, u[i] - (mu+k) + cp[i-1] )
  cm[i] <- max(0, (mu+k) - u[i] + cm[i-1] )
  np[i] <- if_else(cp[i]==0, 0, np[i-1] + 1 )
  nm[i] <- if_else(cm[i]==0, 0, nm[i-1] + 1 )
}

cusum <- tibble(i = 1:(nbre-2), u_i = round(u,2), C_plus = round(cp,2), N_plus = np, C_moins = round(cm,2), N_moins = nm)

kable(cusum)
i u_i C_plus N_plus C_moins N_moins
1 0.42 0.32 1 0.00 0
2 0.45 0.68 2 0.00 0
3 0.01 0.58 3 0.09 1
4 -1.32 0.00 0 1.51 2
5 0.15 0.05 1 1.46 3
6 1.45 1.40 2 0.12 4
7 0.52 1.81 3 0.00 0
8 -1.59 0.12 4 1.69 1
9 -0.77 0.00 0 2.56 2
10 -1.50 0.00 0 4.16 3
11 0.04 0.00 0 4.21 4
12 0.02 0.00 0 4.29 5
13 0.99 0.89 1 3.41 6
14 -0.02 0.77 2 3.53 7
15 1.62 2.29 3 2.00 8
16 0.63 2.82 4 1.48 9
17 -0.15 2.56 5 1.73 10
18 0.24 2.71 6 1.59 11
19 -0.33 2.28 7 2.02 12
20 1.37 3.55 8 0.74 13
21 1.37 4.82 9 0.00 0
22 0.86 5.58 10 0.00 0
23 0.17 5.64 11 0.00 0
24 -1.28 4.26 12 1.38 1
25 1.33 5.49 13 0.15 2
26 0.96 6.36 14 0.00 0
27 -1.40 4.86 15 1.50 1
28 0.70 5.46 16 0.89 2 
plot_cusum <- ggplot( data = cusum %>% select(i, C_plus, C_moins) %>% mutate(C_moins = - C_moins) %>% gather(type, value, C_plus:C_moins) ) +
  aes( x = i, y = value, color = type) +
  geom_line(size = 1) +
  geom_hline(yintercept = 0, color = "black", size = 1) +
  geom_hline(yintercept = h, color = "red", size = 1.2) +
  geom_hline(yintercept = -h, color = "red", size = 1.2)+
  geom_text(x = 1, y = h + 0.5, label = "H", color = "red") +
  geom_text(x = 1, y = -h + 0.5, label = "-H", color = "red") +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:30 ) +
  ylab("C+ / C-") +
  theme_minimal() +
  ggtitle("CUSUM autonome") +
  theme( plot.title = element_text(hjust = 0.5) )

print(plot_cusum)

On constate sur le graphique qu’on dépasse la valeur de H pour la pièce i = 22. En se référant au tableau CUSUM, on lit N+ = 10 lorsque i = 22. On peut donc situer l’évènement (erreur humaine, dérèglement ou autre) qui a provoqué le décalage de la moyenne aux alentours de la pièce i = 12. Cette inférence est satisfaisante, quand on compare à l’origine réelle du problème en i = 12.

4.7 Carte EWMA

Puisque la carte de contrôle EWMA est relativement robuste à la non-normalité, on travaille sur deux échantillons :
- le premier échantillon suit une loi gamma de paramètre shape = 12, scale = 1 - le deuxième échantillon suit une loi gamma de paramètre shape = 15, scale = 1

On génère 20 observations issues du premier échantillon et 20 observations issues du second échantillon. On suppose connues l’espérance et la variance de l’échantillon 1. On construit directement la carte EWMA avec les paramètres lambda = 0.1 et L = 3.

set.seed(3)

nbre_observations <- 40
lambda <- 0.1
l <- 3
shape <- 12
scale <- 1
mu <- shape * scale
shift <- 3
sigma <- sqrt(shape) * scale


big_sample1 <- rgamma( n = 100000, shape = shape, scale = scale )
big_sample2 <- rgamma( n = 100000, shape = shape + shift, scale = scale )
x <- c( sample( big_sample1, nbre_observations/2 ), sample(big_sample2, nbre_observations/2 ) )

z <- vector("double", nbre_observations)
ucl <- vector("double", nbre_observations)
lcl <- vector("double", nbre_observations)

z[1] <- ( lambda * x[1] ) + ( (1-lambda) * mu ) 
for( i in 2:nbre_observations ){
  z[i] <- ( lambda * x[i] ) + ( (1-lambda) * z[i-1] ) 
}

for( i in 1:nbre_observations ){
  ucl[i] <- mu + l*sigma*sqrt( lambda * ( 1 - (1-lambda)**(2*i) ) / (2-lambda) )
  lcl[i] <- mu - l*sigma*sqrt( lambda * ( 1 - (1-lambda)**(2*i) ) / (2-lambda) )
}


ewma_data <- tibble( i = 1:nbre_observations, obs = x, ewma = z, ucl = ucl, lcl = lcl)
ewma_data <- ewma_data %>% gather(data, value, ucl:lcl)

plot_ewma <- ggplot( data = ewma_data ) +
  geom_step( aes( x = i, y = value, color = data) ) +
  geom_line( aes(x = i, y = ewma), color = "black") +
  geom_hline(yintercept = mu, color = "black", size = 1) +
  geom_vline(xintercept = 20.5, linetype = "longdash") +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:nbre_observations ) +
  scale_colour_manual(values = c("red", "blue")) +
  ylab("EWMA") +
  theme_minimal() +
  ggtitle("Carte de contrôle EWMA : loi gamma de paramètres (12,1)",
          subtitle = "lambda = 0.1, L = 3, shift = 3 à partir de l'individu 21.") +
  theme( plot.title = element_text(hjust = 0.5),
         plot.subtitle = element_text(hjust = 0.5) )

print(plot_ewma)

On détecte le shift dans l’espérance à partir de l’individu 32, l’événement ayant réellement eu lieu pour l’individu 21.

5 TD

5.1 Exercice 2

On importe les données, on calcule la moyenne et l’étendue de chaque sous-échantillon de taille 4. Ensuite, on réalise les cartes de contrôle X, R en prenant pour consigne la moyenne des moyennes/étendues de tous les sous-échantillons.

# Import
don <- tibble(x1 = c(6,10,7,8,9,12,16,7,9,15,8,6,16,7,11,15,9,15,8,3),
              x2 = c(9,4,8,9,10,11,10,5,7,16,12,13,9,13,7,10,8,7,6,14),
              x3 = c(10,6,10,6,7,10,8,10,8,10,14,9,13,10,10,11,12,10,9,11),
              x4 = c(15,11,5,13,13,10,9,4,12,13,16,11,15,12,16,14,10,11,12,15)
              ) %>%
  rowwise() %>%
  mutate( xb = mean( c(x1, x2, x3, x4) ), #Moyenne
          rb = max( c(x1, x2, x3, x4)) - min(c(x1, x2, x3, x4) ) ) # Étendue

xbb <- don %>% pull(xb) %>% mean() # Moyenne globale
rbb <- don %>% pull(rb) %>% mean() # Étdendue globale

# Construction des bornes des cartes de contrôle. On lit les coefficients A2, D3 et D4 dans les tables.
a2 <- 0.729
d3 <- 0
d4 <- 2.282 
uclxb <- xbb + a2 * rbb
lclxb <- xbb - a2 * rbb
uclrb <- d4 * rbb
lclrb <- d3 * rbb

# Représentation des cartes de contrôle
plotxb <- ggplot(don) +
  aes( x = 1:nrow(don), y = xb) +
  geom_line(size = 1.2) +
  geom_hline(yintercept = xbb, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclxb, color = "red", size = 1) +
  geom_hline(yintercept = uclxb, color = "red", size = 1)+
  geom_text(x = 2, y = xbb + 0.5, label = "X\u0305\u0305", color = "chartreuse4") +
  geom_text(x = 2, y = lclxb + 0.5, label = "LCL", color = "red") +
  geom_text(x = 2, y = uclxb - 0.5, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = 1:20 ) +
  ylab("X\u0305") +
  theme_minimal() +
  ggtitle("Contrôle de la moyenne X\u0305") +
  theme( plot.title = element_text(hjust = 0.5) )

plotrb <- ggplot(don) +
  aes( x = 1:nrow(don), y = rb) +
  geom_line(size = 1.2) +
  geom_hline(yintercept = rbb, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclrb, color = "red", size = 1) +
  geom_hline(yintercept = uclrb, color = "red", size = 1)+
  geom_text(x = 2, y = rbb + 1, label = "R\u0305\u0305", color = "chartreuse4") +
  geom_text(x = 2, y = lclrb + 1, label = "LCL", color = "red") +
  geom_text(x = 2, y = uclrb - 1, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = 1:20 ) +
  ylab("R\u0305") +
  theme_minimal() +
  ggtitle("Contrôle de l'étendue R") +
  theme( plot.title = element_text(hjust = 0.5) )
  
print(plotxb)

print(plotrb)

On extrait l’intégralité des données dans un seul vecteur afin de réaliser un test de Shapiro-Wilk, et un Q-Q plot.

don_vector <- tibble( x = c( pull(don, x1), pull(don, x2), pull(don, x3), pull(don, x4) ) )

don_vector %>% pull(x) %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.97249, p-value = 0.08182
qqplot <- ggplot(data = don_vector, aes(sample = x)) +
  theme_minimal() +
  stat_qq() +
  stat_qq_line() +
  ggtitle("Q-Q Plot des données")

plot(qqplot)

La p-value est faible mais au-dessus des 5 %. Le test conclut à la normalité de l’échantillon.
Le QQ-plot confirme la normalité de l’échantillon car on peut penser que si les valeurs étaient parfaitement mesurées, alors la distribution serait très proche d’une distribution normale. C’est certainement le manque de précision dans les mesures qui provoquent la présence de nombreux doublons et qui dégradent la normalité de l’échantillon (et donc la p-value du test ci-dessus).

5.2 Exercice 6

On importe les données, on calcule la somme de toutes les erreurs pour chaque jour considéré. Ensuite, on réalise la carte de contrôle c en prenant pour consigne la moyenne des sommes des erreurs.

# Import
don <- tibble(x1 = c(8,11,1,3,3,6,8,4,1,15,1,6,7,2,6,2,11,5,6,2,7,4,4,15,2),
              x2 = c(7,1,1,2,2,3,8,10,6,1,7,7,6,9,14,9,1,5,15,7,5,3,1,2,15),
              x3 = c(1,11,8,5,13,3,2,2,1,3,13,9,3,3,7,4,1,19,5,9,6,8,4,7,3),
              x4 = c(11,2,2,1,6,3,1,6,3,2,5,3,3,8,1,2,3,1,6,2,14,1,20,10,11),
              x5 = c(17,9,5,4,5,1,5,4,2,8,1,1,1,7,8,1,2,3,6,8,10,2,5,17,2)
              ) %>%
  rowwise() %>%
  mutate( cb = sum( c(x1, x2, x3, x4, x5) ) ) # Somme des erreurs journalières



# Construction de la consigne et des LCL/UCL :

cbb <- don %>% pull(cb) %>% mean() # Moyenne de la somme des erreurs journalières
lclc <- cbb - 3 * sqrt( cbb )
uclc <- cbb + 3 * sqrt( cbb )

# Représentation de la carte c

plotc<- ggplot(data = don) +
  aes( x = 1:nrow(don), y = cb) +
  geom_line(size = 1.2) +
  geom_hline(yintercept = cbb, color = "chartreuse4", size = 1.3 ) +
  geom_hline(yintercept = lclc, color = "red", size = 1) +
  geom_hline(yintercept = uclc, color = "red", size = 1)+
  geom_text(x = 2, y = cbb + 2, label = "Consigne", color = "chartreuse4") +
  geom_text(x = 2, y = lclc + 2, label = "LCL", color = "red") +
  geom_text(x = 2, y = uclc + 2, label = "UCL", color = "red") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = seq(0,30,5) ) +
  ylab("c") +
  theme_minimal() +
  ggtitle("Carte c pour le nombre total d'erreurs") +
  theme( plot.title = element_text(hjust = 0.5) )
  
print(plotc)

La carte c produit une alerte le jour 24 : le nombre d’erreurs n’est plus sous contrôle. On ne peut produire de carte t, car pour chaque sous-échantillon on a seulement le nombre total d’erreurs. Il manque le nombre de fiches clients testées pour trouver le nombre d’erreurs. Un exemple de carte t, appelée aussi carte g et h, est proposé dans la section “Cartes de contrôles” ci-dessus.

5.3 Exercice 7

On réalise une carte de contrôle p avec taille d’échantillon variable afind e suivre la qualité du travail du groupe de maintenance. On importe les données, on calcule la moyenne de la fraction de demande de deuxième visite parmi toutes les réparations effectuées. Cette donnée p jouera le rôle de consigne. Puis on calcule les limites LCL/UCL qui varient en fonction de la taille de l’échantillon.

# Import
don <- tibble(n = c(200,250,250,250,200,200,150,150,150,150,100,100,100,200,200,200,200,200,250,250),
              d = c(6,8,9,7,3,4,2,1,0,2,1,0,1,4,5,3,10,4,7,6)
              )

pb <- sum(pull(don, d)) / sum(pull(don,n))
ntot <- nrow(don)

don <- don %>% mutate(
  i = 1:ntot,
  p = d / n,
  sd = sqrt( p * (1-p) / n ),
  lcl = if_else( pb - 3 * sd <0 | d ==0, 0, pb - 3 * sd ),
  ucl = pb + 3 * sd
)

don <- don %>% select(i, p, lcl, ucl) %>% gather(data, value, lcl:ucl)

plotp <- ggplot( data = don ) +
  geom_step( aes( x = i, y = value, color = data) ) +
  geom_line( aes(x = i, y = p), color = "black") +
  geom_hline(yintercept = pb, color = "black", size = 1) +
  geom_vline(xintercept = 20.5, linetype = "longdash") +
  xlab("Échantillon") +
  scale_x_continuous( breaks = 1:ntot ) +
  scale_colour_manual(values = c("red", "blue")) +
  ylab("p") +
  theme_minimal() +
  ggtitle("Carte de contrôle p",
          subtitle = "Taille d'échantillon variable") +
  theme( plot.title = element_text(hjust = 0.5),
         plot.subtitle = element_text(hjust = 0.5) )

print(plotp)

La qualité de la maintenance est globalement sous contrôle.

5.4 Exercice 9

5.5 Exercice 13

don_vector <- tibble(
  x = c(20.0106,20.0090,20.0067,19.9772,20.0001,19.9940,19.9876,20.0042,19.9986,19.9958,20.0075,20.0018,20.0059,19.9975,20.0089,20.0045,19.9891,19.9956,19.9884,20.0154,20.0056,19.9831,20.0040,20.0006,20.0047)
  )

qqplot <- ggplot(data = don_vector, aes(sample = x)) +
  theme_minimal() +
  stat_qq() +
  stat_qq_line() +
  ggtitle("Q-Q Plot des données")

plot(qqplot)

La normalité des données semble être une hypothèse acceptable d’après le Q-Q Plot. Effectuons un test de Shapiro-Wilk :

don_vector %>% pull(x) %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.95029, p-value = 0.2545

Le test accepte l’hypothèse de normalité avec une p-value correcte.
Estimons la moyenne X de l’échantillon,son écart-type s ainsi que les valeurs X + 3s / X - 3s :

xb <- don_vector %>% pull(x) %>% mean()
s <- don_vector %>% pull(x) %>% sd()
l <- xb - 3*s
u <- xb + 3*s

En l’absence de spécifications sur les données, on fait l’estimation de capabilité suivante, basée sur la normalité supposée du processus : 99.73 % des hauteurs de disque se trouveront dans l’intervalle [X - 3s ; X + 3s], i.e. [19.9721, 20.0276].

6 TD 1

6.1 Exercice 1

6.1.1 1) Carte CUSUM

On construit la carte CUSUM associée aux données avec les valeurs classique :
\(k = \frac{1}{2}, K = k \sigma\)
\(h = 4, H = h \sigma\)

don <- tibble( x = c(1045,1055,1037,1064,1095,1008,1050,1087,1125,1146,1139,1169,1151,1128,1238,1125,1163,1188,1146,1167) )
donv <- don %>% pull(x)

mu <- 1050
sigma <- 25
k <- 0.5
h <- 4
K <- k * sigma
H <- h * sigma
nbre <- nrow(don)

cp <- vector("double", nbre)
cm <- vector("double", nbre)
np <- vector("double", nbre)
nm <- vector("double", nbre)

cp[1] <- max(0, donv[1] - (mu+K) )
cm[1] <- max(0, (mu+K) - donv[1] )
np[1] <- if_else(cp[1]==0, 0, 1)
nm[1] <- if_else(cm[1]==0, 0, 1)
  
for(i in 2:nbre){
  cp[i] <- max(0, donv[i] - (mu+K) + cp[i-1] )
  cm[i] <- max(0, (mu+K) - donv[i] + cm[i-1] )
  np[i] <- if_else(cp[i]==0, 0, np[i-1] + 1 )
  nm[i] <- if_else(cm[i]==0, 0, nm[i-1] + 1 )
}

cusum <- tibble(i = 1:nbre, x_i = donv, C_plus = round(cp,2), N_plus = np, C_moins = round(cm,2), N_moins = nm)

kable(cusum)
i x_i C_plus N_plus C_moins N_moins
1 1045 0.0 0 17.5 1
2 1055 0.0 0 25.0 2
3 1037 0.0 0 50.5 3
4 1064 1.5 1 49.0 4
5 1095 34.0 2 16.5 5
6 1008 0.0 0 71.0 6
7 1050 0.0 0 83.5 7
8 1087 24.5 1 59.0 8
9 1125 87.0 2 0.0 0
10 1146 170.5 3 0.0 0
11 1139 247.0 4 0.0 0
12 1169 353.5 5 0.0 0
13 1151 442.0 6 0.0 0
14 1128 507.5 7 0.0 0
15 1238 683.0 8 0.0 0
16 1125 745.5 9 0.0 0
17 1163 846.0 10 0.0 0
18 1188 971.5 11 0.0 0
19 1146 1055.0 12 0.0 0
20 1167 1159.5 13 0.0 0 
plot_cusum <- ggplot( data = cusum %>% select(i, C_plus, C_moins) %>% mutate(C_moins = - C_moins) %>% gather(type, value, C_plus:C_moins) ) +
  aes( x = i, y = value, color = type) +
  geom_line(size = 1) +
  geom_hline(yintercept = 0, color = "black", size = 1) +
  geom_hline(yintercept = H, color = "red", size = 1.2) +
  geom_hline(yintercept = -H, color = "red", size = 1.2)+
  geom_text(x = 1, y = H + 50, label = "H", color = "red") +
  geom_text(x = 1, y = -H + 30, label = "-H", color = "red") +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:nbre ) +
  ylab("C+ / C-") +
  theme_minimal() +
  ggtitle("CUSUM chart") +
  theme( plot.title = element_text(hjust = 0.5) )

print(plot_cusum)

La valeur H est dépassée pour l’observation n°10. On lit N+ = 3 pour i =10, donc on peut situer l’évènement qui a provoqué la perte de contrôle aux alentours de l’observation n°7.

6.1.2 2) Carte EWMA

lambda <- 0.1
l <- 2.7

z <- vector("double", nbre)
ucl <- vector("double", nbre)
lcl <- vector("double", nbre)

z[1] <- ( lambda * donv[1] ) + ( (1-lambda) * mu ) 
for( i in 2:nbre ){
  z[i] <- ( lambda * donv[i] ) + ( (1-lambda) * z[i-1] ) 
}

for( i in 1:nbre ){
  ucl[i] <- mu + l*sigma*sqrt( lambda * ( 1 - (1-lambda)**(2*i) ) / (2-lambda) )
  lcl[i] <- mu - l*sigma*sqrt( lambda * ( 1 - (1-lambda)**(2*i) ) / (2-lambda) )
}


ewma_data <- tibble( i = 1:nbre, obs = donv, ewma = z, ucl = ucl, lcl = lcl)
ewma_data <- ewma_data %>% gather(data, value, ucl:lcl)

plot_ewma <- ggplot( data = ewma_data ) +
  geom_step( aes( x = i, y = value, color = data) ) +
  geom_line( aes(x = i, y = ewma), color = "black") +
  geom_hline(yintercept = mu, color = "black", size = 1) +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:nbre ) +
  scale_colour_manual(values = c("red", "blue")) +
  ylab("EWMA") +
  theme_minimal() +
  ggtitle("Carte de contrôle EWMA",
          subtitle = "lambda = 0.1, L = 2.7") +
  theme( plot.title = element_text(hjust = 0.5),
         plot.subtitle = element_text(hjust = 0.5) )

print(plot_ewma)

Comme la carte CUSUM, la carte EWMA émet aussi une alerte pour la pièce n°10.

6.1.3 3) Carte MA

w <- 6
lcl <- mu - ( 3 * sigma / sqrt(w) )
ucl <- mu + ( 3 * sigma / sqrt(w) )

m <- vector("double", nbre)
m[1] <- donv[1]
for(i in 2:nbre){
  if(i <= w){
    m[i] <- mean(donv[1:i])
  }
  else{
    m[i] <- mean(donv[(i-w+1):i])
  }
}

l <- vector("double", nbre)
u <- vector("double", nbre)
for(i in 1:nbre){
  if(i<w){
    l[i] <- mu - ( 3 * sigma / sqrt(i) )
    u[i] <- mu + ( 3 * sigma / sqrt(i) )
  }else{
    l[i] <- lcl
    u[i] <- ucl
  }
}

ma_data <- tibble( i = 1:nbre, obs = donv, ma = m, ucl = u, lcl = l)
ma_data <- ma_data %>% gather(data, value, ucl:lcl)

plot_ma <- ggplot( data = ma_data ) +
  geom_step( aes( x = i, y = value, color = data) ) +
  geom_line( aes(x = i, y = ma), color = "black") +
  geom_hline(yintercept = mu, color = "black", size = 1) +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:nbre ) +
  scale_colour_manual(values = c("red", "blue")) +
  ylab("MA") +
  theme_minimal() +
  ggtitle("Carte de contrôle MA",
          subtitle = "w = 6") +
  theme( plot.title = element_text(hjust = 0.5),
         plot.subtitle = element_text(hjust = 0.5) )

print(plot_ma)

La carte MA, comme les cartes précédentes, émet une alerte pour l’observation n°10.

6.2 Exercice 2

6.2.1 1) Carte CUSUM

On construit la carte CUSUM associée aux données et aux valeurs spécifiées :

don <- tibble( x = c(8.00,8.01,8.02,8.01,8.00,8.01,8.06,8.07,8.01,8.04,8.02,8.01,8.05,8.04,8.03,8.05,8.06,8.04,8.05,8.06,8.04,8.02,8.03,8.05) )
donv <- don %>% pull(x)

mu <- 8.02
sigma <- 0.05
k <- 0.5
h <- 4.77
K <- k * sigma
H <- h * sigma
nbre <- nrow(don)

cp <- vector("double", nbre)
cm <- vector("double", nbre)
np <- vector("double", nbre)
nm <- vector("double", nbre)

cp[1] <- max(0, donv[1] - (mu+K) )
cm[1] <- max(0, (mu+K) - donv[1] )
np[1] <- if_else(cp[1]==0, 0, 1)
nm[1] <- if_else(cm[1]==0, 0, 1)
  
for(i in 2:nbre){
  cp[i] <- max(0, donv[i] - (mu+K) + cp[i-1] )
  cm[i] <- max(0, (mu+K) - donv[i] + cm[i-1] )
  np[i] <- if_else(cp[i]==0, 0, np[i-1] + 1 )
  nm[i] <- if_else(cm[i]==0, 0, nm[i-1] + 1 )
}

cusum <- tibble(i = 1:nbre, x_i = donv, C_plus = round(cp,2), N_plus = np, C_moins = round(cm,2), N_moins = nm)

kable(cusum)
i x_i C_plus N_plus C_moins N_moins
1 8.00 0.00 0 0.04 1
2 8.01 0.00 0 0.08 2
3 8.02 0.00 0 0.11 3
4 8.01 0.00 0 0.14 4
5 8.00 0.00 0 0.19 5
6 8.01 0.00 0 0.22 6
7 8.06 0.02 1 0.21 7
8 8.07 0.04 2 0.18 8
9 8.01 0.01 3 0.21 9
10 8.04 0.00 0 0.22 10
11 8.02 0.00 0 0.25 11
12 8.01 0.00 0 0.28 12
13 8.05 0.01 1 0.28 13
14 8.04 0.00 0 0.28 14
15 8.03 0.00 0 0.30 15
16 8.05 0.01 1 0.29 16
17 8.06 0.02 2 0.28 17
18 8.04 0.02 3 0.28 18
19 8.05 0.02 4 0.28 19
20 8.06 0.04 5 0.26 20
21 8.04 0.03 6 0.27 21
22 8.02 0.01 7 0.29 22
23 8.03 0.00 0 0.31 23
24 8.05 0.01 1 0.30 24 
plot_cusum <- ggplot( data = cusum %>% select(i, C_plus, C_moins) %>% mutate(C_moins = - C_moins) %>% gather(type, value, C_plus:C_moins) ) +
  aes( x = i, y = value, color = type) +
  geom_line(size = 1) +
  geom_hline(yintercept = 0, color = "black", size = 1) +
  geom_hline(yintercept = H, color = "red", size = 1.2) +
  geom_hline(yintercept = -H, color = "red", size = 1.2)+
  geom_text(x = 1, y = H + 1, label = "H", color = "red") +
  geom_text(x = 1, y = -H + 1, label = "-H", color = "red") +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:nbre ) +
  ylab("C+ / C-") +
  theme_minimal() +
  ggtitle("CUSUM chart") +
  theme( plot.title = element_text(hjust = 0.5) )

print(plot_cusum)

La valeur observée pour i = 11 est inférieure à -H : la carte émet une alerte. On lit N- = 11 pour i = 11, donc on peut penser que l’évènement qui a provoqué la perte de contrôle s’est produit avant même le début des mesures considérées ici.

6.2.2 2) Test du \(\chi^{2}\) d’égalité de \(\sigma\) à une valeur donnée (0.05) :

Ce test suppose que les données suivent une loi normale, testons d’abord cela :

shapiro.test(donv)
## 
##  Shapiro-Wilk normality test
## 
## data:  donv
## W = 0.93407, p-value = 0.1202

Le test de Shapiro-Wilk accepte l’hypothèse de normalité.

sigma2 <- sigma**2
varTest(donv, alternative = "two.sided", conf.level = 0.95, sigma.squared = sigma2)
## 
##  Chi-Squared Test on Variance
## 
## data:  donv
## Chi-Squared = 4.1, df = 23, p-value = 8.636e-06
## alternative hypothesis: true variance is not equal to 0.0025
## 95 percent confidence interval:
##  0.0002692011 0.0008769264
## sample estimates:
##     variance 
## 0.0004456522

L’égalité \(\sigma = 0.05\) est rejetée avec une p-value très faible.

6.2.3 3) Carte CUSUM : 2ème épisode

On construit la carte CUSUM associée aux données et aux valeurs spécifiées :

mu <- 8.01
k <- 0.25
K <- k * sigma

cp <- vector("double", nbre)
cm <- vector("double", nbre)
np <- vector("double", nbre)
nm <- vector("double", nbre)

cp[1] <- max(0, donv[1] - (mu+K) )
cm[1] <- max(0, (mu+K) - donv[1] )
np[1] <- if_else(cp[1]==0, 0, 1)
nm[1] <- if_else(cm[1]==0, 0, 1)
  
for(i in 2:nbre){
  cp[i] <- max(0, donv[i] - (mu+K) + cp[i-1] )
  cm[i] <- max(0, (mu+K) - donv[i] + cm[i-1] )
  np[i] <- if_else(cp[i]==0, 0, np[i-1] + 1 )
  nm[i] <- if_else(cm[i]==0, 0, nm[i-1] + 1 )
}

cusum <- tibble(i = 1:nbre, x_i = donv, C_plus = round(cp,2), N_plus = np, C_moins = round(cm,2), N_moins = nm)

kable(cusum)
i x_i C_plus N_plus C_moins N_moins
1 8.00 0.00 0 0.02 1
2 8.01 0.00 0 0.03 2
3 8.02 0.00 0 0.04 3
4 8.01 0.00 0 0.05 4
5 8.00 0.00 0 0.07 5
6 8.01 0.00 0 0.08 6
7 8.06 0.04 1 0.05 7
8 8.07 0.09 2 0.00 0
9 8.01 0.07 3 0.01 1
10 8.04 0.09 4 0.00 0
11 8.02 0.09 5 0.00 1
12 8.01 0.08 6 0.01 2
13 8.05 0.10 7 0.00 0
14 8.04 0.12 8 0.00 0
15 8.03 0.13 9 0.00 0
16 8.05 0.16 10 0.00 0
17 8.06 0.19 11 0.00 0
18 8.04 0.21 12 0.00 0
19 8.05 0.24 13 0.00 0
20 8.06 0.28 14 0.00 0
21 8.04 0.29 15 0.00 0
22 8.02 0.29 16 0.00 1
23 8.03 0.30 17 0.00 0
24 8.05 0.33 18 0.00 0 
plot_cusum <- ggplot( data = cusum %>% select(i, C_plus, C_moins) %>% mutate(C_moins = - C_moins) %>% gather(type, value, C_plus:C_moins) ) +
  aes( x = i, y = value, color = type) +
  geom_line(size = 1) +
  geom_hline(yintercept = 0, color = "black", size = 1) +
  geom_hline(yintercept = H, color = "red", size = 1.2) +
  geom_hline(yintercept = -H, color = "red", size = 1.2)+
  geom_text(x = 1, y = H + 1, label = "H", color = "red") +
  geom_text(x = 1, y = -H + 1, label = "-H", color = "red") +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:nbre ) +
  ylab("C+ / C-") +
  theme_minimal() +
  ggtitle("CUSUM chart") +
  theme( plot.title = element_text(hjust = 0.5) )

print(plot_cusum)

6.2.4 4) et 5) Comparaison des résultats

La première carte émet une alerte pour i = 11 car C- < -H. La deuxième carte émet une alerte pour i = 19 car C+ > H.
On peut penser que la première carte n’est pas fiable du tout car \(\sigma\) est sur-estimé (en réalité \(\sigma < 0.05\)). Il en résulte que le paramètre K est trop grand, et que la quantité \((\mu_{0}-K)-x_{i}\) est presque toujours négative, et qu’on a une abondance anormale de C- non-nuls. La deuxième carte paraît plus fiable car on a diminué le pramètre k (k=0.25), ce qui va dans le sens d’une variance plus faible et donc d’une carte plus précise.

6.2.5 6) Carte EWMA

lambda <- 0.2
l <- 3

z <- vector("double", nbre)
ucl <- vector("double", nbre)
lcl <- vector("double", nbre)

z[1] <- ( lambda * donv[1] ) + ( (1-lambda) * mu ) 
for( i in 2:nbre ){
  z[i] <- ( lambda * donv[i] ) + ( (1-lambda) * z[i-1] ) 
}

for( i in 1:nbre ){
  ucl[i] <- mu + l*sigma*sqrt( lambda * ( 1 - (1-lambda)**(2*i) ) / (2-lambda) )
  lcl[i] <- mu - l*sigma*sqrt( lambda * ( 1 - (1-lambda)**(2*i) ) / (2-lambda) )
}


ewma_data <- tibble( i = 1:nbre, obs = donv, ewma = z, ucl = ucl, lcl = lcl)
ewma_data <- ewma_data %>% gather(data, value, ucl:lcl)

plot_ewma <- ggplot( data = ewma_data ) +
  geom_step( aes( x = i, y = value, color = data) ) +
  geom_line( aes(x = i, y = ewma), color = "black") +
  geom_hline(yintercept = mu, color = "black", size = 1) +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:nbre ) +
  scale_colour_manual(values = c("red", "blue")) +
  ylab("EWMA") +
  theme_minimal() +
  ggtitle("Carte de contrôle EWMA",
          subtitle = "lambda = 0.2, L = 3") +
  theme( plot.title = element_text(hjust = 0.5),
         plot.subtitle = element_text(hjust = 0.5) )

print(plot_ewma)

6.2.6 7) Carte MA

w <- 5
lcl <- mu - ( 3 * sigma / sqrt(w) )
ucl <- mu + ( 3 * sigma / sqrt(w) )

m <- vector("double", nbre)
m[1] <- donv[1]
for(i in 2:nbre){
  if(i <= w){
    m[i] <- mean(donv[1:i])
  }
  else{
    m[i] <- mean(donv[(i-w+1):i])
  }
}

l <- vector("double", nbre)
u <- vector("double", nbre)
for(i in 1:nbre){
  if(i<w){
    l[i] <- mu - ( 3 * sigma / sqrt(i) )
    u[i] <- mu + ( 3 * sigma / sqrt(i) )
  }else{
    l[i] <- lcl
    u[i] <- ucl
  }
}

ma_data <- tibble( i = 1:nbre, obs = donv, ma = m, ucl = u, lcl = l)
ma_data <- ma_data %>% gather(data, value, ucl:lcl)

plot_ma <- ggplot( data = ma_data ) +
  geom_step( aes( x = i, y = value, color = data) ) +
  geom_line( aes(x = i, y = ma), color = "black") +
  geom_hline(yintercept = mu, color = "black", size = 1) +
  xlab("Pièce") +
  scale_x_continuous( breaks = 1:nbre ) +
  scale_colour_manual(values = c("red", "blue")) +
  ylab("MA") +
  theme_minimal() +
  ggtitle("Carte de contrôle MA",
          subtitle = "w = 5") +
  theme( plot.title = element_text(hjust = 0.5),
         plot.subtitle = element_text(hjust = 0.5) )

print(plot_ma)

Les cartes EWMA et MA ne produisent pas d’alerte, mais les graphiques indiquent malgré tout une tendance à la hausse pour la moyenne. Cela confirme que la carte CUSUM n°2 (question 3) était plus adapté que la carte CUSUM n°1 (question 1.).

6.3 Exercice 3

7 TD 2

7.1 Exercice 1

7.2 Exercice 2

rc <- function(l){
  return(
    ppois(3,l) + sum( ppois(5:0,l) * (ppois(4:9,l)-ppois(3:8,l))  )
    )
}

donc <- tibble( x = seq(8,20,0.01)) %>% rowwise() %>% mutate( p = rc(x) )

risque_client <- ggplot(data = donc) +
  aes( x = x, y = p ) +
  geom_line( size = 1.1, color = "black" ) +
  geom_point( x = 8, y = rc(8), color = 'red', size = 3 ) +
  scale_x_continuous(breaks = 8:20) +
  xlab( expression(lambda) ) +
  ylab("Risque client") +
  ggtitle("Risque client avec approximation de Poisson") +
  theme_minimal()

print(risque_client)

print(rc(8))
## [1] 0.07001466
rf <- function(l){
  return(
    1 - ppois(9,l) + sum( ( 1 - ppois(5:0,l) ) * (ppois(4:9,l)-ppois(3:8,l))  )
  )
}

donf <- tibble( x = seq(0,3,0.01)) %>% rowwise() %>% mutate( p = rf(x) )

risque_fournisseur <- ggplot(data = donf) +
  aes( x = x, y = p ) +
  geom_line( size = 1.1, color = "black" ) +
  geom_point( x = 3, y = rf(3), color = 'red', size = 3 ) +
  scale_x_continuous(breaks = 0:3) +
  xlab( expression(lambda) ) +
  ylab("Risque fournisseur") +
  ggtitle("Risque fournisseur avec approximation de Poisson") +
  theme_minimal()

print(risque_fournisseur)

print(rf(3))
## [1] 0.07312669
n3 <- function(l){
  return(
    100 * ( ppois(3,l) + 1 - ppois(9,l) ) + 200 * ( ppois(9,l) - ppois(4,l) )
  )
}

donn <- tibble( x = seq(0,20,0.01)) %>% rowwise() %>% mutate( p = n3(x) )

esperance_n2 <- ggplot(data = donn) +
  aes( x = x, y = p ) +
  geom_line( size = 1.1, color = "black" ) +
  #geom_point( x = 8, y = rc(8), color = 'red', size = 3 ) +
  scale_x_continuous(breaks = 0:20) +
  xlab( expression(lambda) ) +
  ylab("Espérance du nombre de pièces prélevées") +
  ggtitle("Espérance du nombre de pièces prélevées dans le cas n°3") +
  theme_minimal()

print(esperance_n2)

7.3 Exercice 4